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Colloidal suspensions have proven lo be excellent 
model systems for the study of condensed matter and its 
phase behavior [11. Many of the properties of colloidal 
suspensions can be investigated with a systematic varia- 
tion of the characteristics of the systems and, in addition, 
the energy, length and lime scales associated with them 
allow for experimental probing of otherwise inacces- 
sible regimes [2, 3]. The latter property also makes 
colloidal systems vulnerable to external influences such 
as gravity. Experiments performed in micro-gravity by 
Chaikin and Russell have been invaluable in extracting 
the true behavior of the systems without an external 
field [3]. Wcitz and Pusey (private communication) in- 
tend to use mixtures of colloidal particles with additives 
such as polymers to induce aggregation and form weak, 
tenuous, highly disordered fractal structures that would 
be stable in the absence of gravitational forces. 

Between any two colloidal particles, there is an at- 
tractive force caused by the interactions between the 
fluctuating dipole moments of their constituent atoms. 
This van der Waals force leads to irreversible aggre- 
gation of unprotected particles and it is desirable to 
diminish its magnitude [4] and provide a stabilization 
mechanism. Common mechanisms are steric stabiliza- 
tion and charge stabilization [1]. For steric stabilization, 
colloidal particles are coated with polymers so that 
there is a short-range highly repulsive (nearly hard core) 
force between two colloidal particles when their respec- 
tive polymer layers are compressed. We shall consider 
these types of colloidal particles later. 

Charge stabilization occurs when charged surface 
groups on colloidal particles dissociate, ionizing the 
particle and emitting counterions into the solution. The 
standard theory of interaction in these charged colloidal 
suspensions is due to Derjaguin, Landau, Verwey and 
Overbeek (DLVO) [5] and the DLVO potential includes 
a non-negligible screened Coulomb repulsion 


in addition to the short-range van der Waals attraction. 
The inverse Debye screening length, n and the effective 
charge, Z* are given by n 2 = Anpe 2 /ksTe, Z * = 
Z exp(na ) /(l + Ka). Here, e is the dielectric constant of 


the solvent, p is the density of counterions, a is the radius 
of a colloidal “macroion”, Z is its bare charge, k B and 
T are the Boltzmann constant and temperature. The 
screened Coulomb potential is obtained in the mean- 
field linearized Poisson-Boltzmann approximation and 
is thus expected lo have have limited applicability, 
although the functional form of the potential with a 
renormalized charge has been argued to remain valid in 
a wider regime [6]. 

Indeed, the DLVO theory has been shown to explain 
some experimental data [7] but fails to account for others 
[8]. The main experimental challenges to the theory are 
numerous indications of a long-range attraction between 
colloidal particles. The latter are evident in experiments 
where particles arc confined near a wall or in suspensions 
at higher densities. When two isolated particles at a low 
density have been considered, no attraction has been 
found. Hence, many-body effects might be thought to 
be responsible for this unexpected attractive force. 

To study the role of confinement and evaluate many- 
body forces, we consider a system consisting of spherical 
colloidal macroions, a counterion density distribution 
p(f) and a homogeneous solvent of a specified dielectric 
constant. The effective Hamiltonian for our system can 
be expressed as: 

u = \Y J m a 2 + Y u{ ^ j - + T (2) 

/ T,J<J 

Ah and Rj are the masses and coordinates of the 
macroions and U(\Rj - Rj |) is the macroion-macroion 
interaction energy. In the expression for the free energy 
of the counterions [9, 10], 

T — J d*rp(r)[(p eT t + — 0uif ^ "f < ^ corr (p)] (^) 

we include: the macroion-counterion and counterion- 
wall interaction potentials, 4> ex t\ the Coulomb interac- 
tions between counterions, <f>i U i (r) — / dV ; ideal 
gas, fad, and correlation, contributions. The ideal 
gas term is <j> l( t = kBT\og[\%p(r)] - 1, where A b is the 
de Broglie thermal wavelength and we employ the free 
energy of a one-component plasma [11] in evaluating 
<t> C orr . We use pseudopotentials for macroion-counterion 
interactions [10] to suppress the large variation of the 
counterion density near a macroion core. 
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Given T, the forces between macroions due to coun- 
terions can be evaluated using the Hellmann-Feynman 
theorem. The total force acting on a macroion becomes: 

f, = -v s T( P (7),<l>e X t) -v ~ R>\)- W 

J 

The equilibrium counterion density, p(f ), itself can be 
obtained from a functional minimization: ^ = 0 . An 
efficient scheme for solving a functional minimization 
of this type was developed by Car and Parrinello [12] in 
the context of quantum-mechanical electronic properties 
calculations and also used by Lowen, Madden and 
Hansen [10] for colloidal suspensions. Expressing the 
counterion density in terms of a wave function y j (f), 

p(f) = |y>(f)| 2 , one can consider the Fourier components 

of v>(r), i/’£, to be dynamical variables in a Lagrangian: 

£ - \ £ mg I4l 2 +s E A l,R, 2 -r- E U(\Rj-R,\) 

(5) 

with a constraint: f d 3 rp(r) = NZ. N is the number 
of macroions, m % are Fictitious masses determining the 
time-scale of the dynamics of counterion density. We 
used this Lagrangian to solve the functional minimiza- 
tion by performing dynamical simulated annealing on 
ihj. to obtain the minimum of T for a Fixed macroion 
configuration. 

We started out with a conFiguration of two macroions 
in a periodic cubic box of length l//m at a temperature 
T = 300 1\ in water, e = 78. The counterion density was 
evaluated on a 96 5 grid and the radius and charge of a 
macroion were a = lOOnm and Z = 200c” (e stands for 
the elementary charge). The forces between macroions 
at a scries of different separations were measured. At 
this low packing fraction, relatively low surface charge 
and counterion concentration, the DLVO pair-potential 
between two macroions is expected to hold. We plot 
the DLVO prediction for the periodic system in Fig. 1. 
We show our density functional minimization results in 
the same figure as well as the force resulting from an 
optimal Fitting to a pair-potential of the DLVO form. 
The numerical values for the inverse screening length 
and effective charge for the optimal fit are presented as 
the First set in Table 1 . 

In order to check the validity of pair- wise additivity 
that we had tacitly assumed above, we placed three 
colloidal particles in a unit cell with periodic boundary 
conditions. Maintaining our previous pseudopotential 
parameters, we studied the following geometries. First, 
we positioned the three macroions in configurations of 
equilateral triangles of different sizes and calculated 



Figure 1: The total force on a macroion in a periodic 
system with two macroions and free counterions in the 
primary simulation cell. The DLVO prediction with no 
adjustable parameters (solid line) is compared to nu- 
merical local density approximation data (open circles) 
and the best fit of a DLVO-tvpe potential is also plotted 
(dashed line). The force is scaled to (e“ /nrn) 2 . 

the force on one of the macroions as a function of 
the edge of the triangle. In Fig. 2, we plot our data 
for the force in these conFigurations as a function of 
the distance and compare it to the force obtainable 
from an additive pair-potential. Secondly, we also 
placed the three macroions in an asymmetric triangular 
conFiguration where two macroions were close to each 
other (separation of 20077m) and a third macroion was at 
a distance of 473nm from either of the two. We measure 
a force of 1.24 - 10 -2 (e“ Inin) 2 ±0.05 * 10 “ (c /nm)~ on 
either one of the First macroions while pairwise additive 
forces predict 1 .22- 10 -2 (e~ /nm) 2 ±0.05 - 10 2 (e jnm ) 2 
for the value. Thus, we detect no three-body component 
in the effective forces. 

To study forces between macroions in a conFined 
system, we introduced two parallel short-range repulsive 
walls. We placed the macroions at different separations 
and distances from a wall and measured the effective 
forces between them. We Fit our data to a DLVO-type 
force and the effective parameters obtainable from the 
Fit are summarized in Table 1. Our results can be 
understood as follows. The introduction of repulsive 
walls increased the counterion density in the middle 
of the simulation cell, while decreasing the Debye 
screening length in that region (compare sets 1 and 
2 in the table). On the other hand, near the walls, 
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R x {nm) 

KG 

2* (e") 

Z (e-) 

1 

- 

0.374 

212.5 

200.9 

2 

500.0 

0.454 

220.2 

203.3 

3 

350.0 

0.444 

220.5 

204.2 

4 

200.0 

0.308 

213.8 

205.5 


Table 1: The parameters corresponding to the optimal 
effective pair interaction potential between macroions 
in different geometries. The first set of data was ob- 
tained with periodic boundary conditions while walls 
were introduced in the last three. R x measures the dis- 
tance from a wall, kci and Z * are obtained from the fit 5 
the bare charge Z is calculated from na and Z * . The 
DLVO prediction for a periodic system is: Z = 200, 
Z‘ = 205.8 and na = 0.258. 

the counterion density must decrease and the Debye 
screening length consequently increase — the latter can 
be seen by comparing sets 2, 3 and 4. Thus, this simple 
model for confinement changes the effective parameters 
in the DLVO potential but does not introduce any 
fundamentally new behavior. 

Besides electrostatic interactions, entropic depletion 
effects that arise from (hard-core) exclusion play an 
important role in determining the behavior of multicom- 
ponent colloidal suspensions [13, 14]. It is well known 
that the addition of free polymers to a suspension of 
colloidal panicles can cause an effective attraction. The 
first successful attempt to predict and explain the phe- 
nomenon was made by Asakura and Oosawa (AO) [15]. 
In their approach, there is an exclusion volume around 
each large sphere within which the centers of small 
ones cannot penetrate. As the entirely entropic free 
energy of small spheres depends on the volume acces- 
sible to them, the favored configurations are the ones 
where the exclusion volumes of large spheres overlap. 
This implies that there ought to be a force, pushing the 
large spheres towards each other in order to increase 
the entropy of the small ones. A similar mechanism 
explains the attraction of a large sphere towards a wall. 
Quantitative AO calculations ignore the interactions be- 
tween the small spheres completely and even though 
an overall acceptable experimental agreement with the 
AO prediction has been reported [13, 16], the ideal 
gas approximation based AO theory does not provide a 
satisfactory description of entropic potentials and forces 
in dense colloidal suspensions. 

To go beyond the simple geometric arguments of 
AO requires a detailed theory for the structure of a binary 



Figure 2: Test of the additivity of the effective macroion 
potential. Three macroions were placed in equilateral 
triangular configurations and the force on a macroion 
as a function of the edge of the triangle was measured. 
Open circles are our data whereas the solid line indi- 
cates the force from an additive DLVO- type potential 
with the parameters: Z * = 212. 5c _ , na — 0.374. The 
force is scaled to (f'/nm) 2 . 

fluid. Even though significant progress has been made 
[17], this is a difficult task for an arbitrary binary system. 
Wc consider instead the simplest situation in which all 
the spheres are of the same size and study depletion 
forces and potentials by an analysis of the data lor this 
system. All our numerical data were obtained from 
hard-sphere molecular dynamics simulations where the 
system evolves on a collision-to-collision basis: all 
particles move freely until two of them come in contact, 
then an elastic collision between the two occurs, after 
which all particles move freely until the next collision. 
Desired quantities are measured as equilibrium time 
averages in these simulations. 

Let us first consider the depletion potential between 
two hard spheres in a “solution” of other hard spheres 
of the same size. The potential of mean force, IT , is 
related to the pair correlation function, g(x ), as [ 1 8] 

g(x) = exp(— W/kuT). (6) 

The depletion force therefore turns out to be: 

F/tBT= 2^m, ( 7 ) 

where x is the distance between two particles. An ex- 
ample of the depletion force between two hard spheres 
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Figure 3: The depletion force as a function of the dis- 
tance between two hard spheres in a mixture of other 
hard spheres of the same size at a packing fraction 
0 = 0.45. The distance is measured in hard sphere 
diameters, <7, and the force is scaled to k B Tjc r. 


in a hard-sphere fluid derived from a pair-correlation 
function using Eq. (7) is given in Figure 3. At a packing 
fraction of 4> = 0.45, it suffices to consider 108 hard 
spheres in a cubic box with periodic boundary condi- 
tions to obtain the first features in the force profile. We 
can prove the validity of our calculation by measuring 
the depiction force directly by recording the total mo- 
mentum transfer in a unit time on a sphere at a known 
distance from another one. We compare both of those 
force measurements to the AO prediction in Figure 3. 

Similarly, the depletion force between a sphere and 
a wall can be obtained from the density profile: 


F/ksT 


0 \u(p(x)) 
dx 


(8) 


where r is the direction perpendicular to the wall and 
p(x) is the local density of the hard sphere fluid. The 
well-known layering of fluid molecules near a wall 
leads again to a more complex structure for the force 
than AO predicts. See Fig. 4 for the force data and 
the AO prediction. We also plot the force obtained 
from momentum transfer measurements. The data was 
obtained with the same number of spheres at the same 
density as before. 

Depletion forces in complex geometries with step 
edges and non-trivial curvature have been measured and 
argued to have profound consequences for cellular biol- 
ogy or entropic control and directed motion of colloidal 
particles. The AO-type calculations of entropic forces 


Figure 4: The depletion force between a hard sphere 
and a wall for the system described in the text, x is mea- 
sured in hard sphere diameters and the force is scaled to 
k B T/a . The walls are parallel to the yz plane and are 
located at x — i2.37. In the simulation, presented in 
this figure, thermal walls were used but reflecting walls 
result in a qualitatively similar force profile. 


and potentials in these geometries arc simple enough 
to carry out while more precise predictions quickly 
become rather laborious and, to our knowledge, have 
not been attempted. To go beyond the AO theory, we 
calculated entropic forces and potentials for a system of 
324 hard spheres at a packing fraction of about 0.45 in 
a three-dimensional T-shaped channel. Potentials near 
step-edges and comers are of principal interest and are 
shown in Figure 5. A comparison with the AO predic- 
tion is also shown in the Figure. The potential barrier 
repelling a hard sphere from a step-edge is simply a re- 
flection of the density decrease near it. Analogously, the 
potential minimum attracting particles towards a comer 
and the sharp increase in density there are different 
manifestations of the same effect. 

There is an interesting and important feature to the 
depletion forces that already manifests itself in the AO 
approach. Depletion forces need not be pairwise ad- 
ditive. In the AO approximation, the force between 
two particles arises from the overlap of excluded vol- 
umes. Three-body effects become relevant when the 
hard spheres are considered in configurations in which 
their pairwise overlap volumes themselves overlap, re- 
sulting in an over-counting of the actual overlap volume. 
An example could be a triangular configuration. The AO 
prediction for an equilateral configuration is calculable. 
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Figure 5: The effective potential profiles for a hard 

sphere fluid confined in a T-shaped channel. To form 
the channel, walls in the yz plane were placed at: 
x = 4.74a; x = 0 for |y| > 2.37a; x = -4.74a for 
\y\ < 2.37a; and in the xz plane at y — ±2. 37a for 
—4.74a < x < 0. (Periodic boundary conditions were 
used otherwise.) The upper figure displays the effective 
potential at x — 0, i. e. along a “step,” the edges of 
which are located at y = ±2.37a. The lower figure is 
the effective potential at x = —4.74a where the “cor- 
ners” axe located at y = ±2. 37a. The potential is given 
in units of thermal energy and y is measured in hard 
sphere diameters. 


For the packing fraction 0 = 0.45, the AO two-body 
force at two sphere contact is Fao = —2.0 k B T/a. If 
forces were additive, the magnitude of the total force 
on a sphere in this triangular configuration would be 
2F ao cos( 7 t/ 6 ) = -3.5 k B T/a. Instead, when all over- 
lap volumes arc correctly accounted for, an AO- like 
analysis predicts instead a force of —2.9 k B T/<r. The 
exact three-body force can in principle be found from 
a three-point correlation function but obtaining reliable 
statistics for this measurement is rather time consuming. 
Instead, we prepared the described triangular configu- 
ration and measured the momentum transfer on each 
of the spheres in a simulation. The two-body force 
al contact for <,6 = 0.45 is F sim = -C.8 k B T/a (see 
Fig. 3). The actual measured force, -10.0 k B T/a, 
is again significantly lower than the additive pairwise 
force: 2F & i m cos(7r/6) = —11.8 k B T/a. Indeed, the 
three-body component in this geometry for our system 
of hard spheres is a significant percentage of the total 
force. 


In conclusion, we have performed calculations of 
effective forces between charged colloidal particles in 
different configurations and studied effective depletion 
forces in hard sphere colloids. In the first case, our 
data are well fit to screened Coulomb type potentials. It 
should be noted that in our calculations, the walls were 
not charged and dielectric discontinuities between the 
suspension and the walls were not considered. One may 
speculate that long-range attraction cannot be accounted 
for by static calculations and that they probably originate 
from temporal fluctuations. In the second case, a com- 
parison of our numerical data and the Asakura-Oosawa 
theory indicates that the AO theory underestimates the 
magnitude of depletion forces al contact, ignores the 
complex nature of entropic interactions, consisting of 
several maxima and minima and underestimates the 
range of the interactions. 
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